Growing spatial correlations of particle displacements in a simulated liquid on cooling 

toward the glass transition 
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We define a correlation function that quantifies the spatial correlation of single-particle displace- 
ments in liquids and amorphous materials. We show for an equilibrium liquid that this function 
is related to fluctuations in a bulk dynamical variable. We evaluate this function using computer 
simulations of an equilibrium glass- forming liquid, and show that long range spatial correlations of 
displacements emerge and grow on cooling toward the mode coupling critical temperature. 



Liquids cooled toward their glass transition exhibit re- 
markable dynamical behavior |lj]. The initial slowing of 
transport processes for liquids at temperatures T well 
above their glass transition temperature T g is described 
to great extent by the mode coupling theory (MCT) 
which predicts diverging relaxation times at the 
mode coupling dynamical critical temperature T c . The 
dynamical singularity of MCT occurs without a diverg- 
ing or even growing static correlation length j^]. Yet 
recent studies show that in the range of T well described 
by MCT, simulated glassforming liquids exhibit spatially 
heterogeneous dynamics p-Q. In this Letter, we define 
a correlation function that quantifies the spatial corre- 
lation of particle displacements and evaluate this func- 
tion for a simulated Lennard-Jones liquid. We find that 
spatial correlations of displacement arise and become in- 
creasingly long ranged on cooling toward T c . 

First, we briefly review the conventional static corre- 
lation function that describes the average microscopic 
structure of a liquid. We use a definition that will 
readily facilitate an extension to a new correlation func- 
tion for particle displacements. Consider a liquid in 
the grand canonical ensemble confined to a volume V, 
consisting of identical particles, each with no internal 
degrees of freedom. Let the position of each parti- 
cle i be denoted r^. In equilibrium the structure of 
a homogeneous liquid can be quantified by G(r), the 
"density-density" correlation function (l^Jl^], defined as 

G(r) = Jdr'/[n(r' + r) - (n)] [n(r') - (n)]^. Here, 

n(r) = Yli=i — r 0' an d (• • •) indicates an ensemble 
average. N — J dr n(r) is the number of particles present 
in a particular configuration. For a homogeneous liquid 
the density p — (n) — (N) /V . If the liquid is isotropic, 
G(r) further reduces simply to G(r), where r = |r|. G(r) 
measures the spatial correlations of fluctuations of local 
density away from the average value. The pair correlation 
function g(r) conventionally presented to characterize the 
structure of a liquid is proportional to the "distinct" part 
of G(r): G(r) = (N)S(r) + (n)(N) [g{r) - 1], where g(r) 
can be written, 
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The Fourier transform of G(r) gives the static structure 
factor 5(q) = ^N^ 1 J2? =1 J2f=i exp[-«q ■ (r, -r^)]^. 

To determine the behavior of G(r) for large r, it is 
useful to evaluate the fluctuations of N, which are re- 
lated to the volume integral of G(r), and also related to 
a thermodynamic response function, the isothermal com- 
pressibility k 



([N- (A0f)= / drG(r) = {n)(N)kTn, 
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where k is Boltzmann's constant. The convergence or di- 
vergence of the volume integral of G(r) depends on how 
rapidly G(r) decays to zero as r — > oo. If the integral 
converges, G(r) is "short ranged" ; if it diverges, G(r) is 
"long ranged". Near a conventional critical point, k di- 
verges, macroscopic density fluctuations occur, and the 
behavior of G(r) approaches that of a long ranged func- 
tion. 

To develop a simple spatial correlation function for a 
local dynamical property in a liquid, we consider for a 
particle i its displacement Hi(t, At) = |rj(i + At) — 
over some interval of time At, starting from a reference 
time t. We examine the spatial correlations of these dis- 
placements by modifying the definition of G(r) so that 
the contribution of a particle i to the correlation function 
is weighted by That is, we define a "displacement- 
displacement" correlation function ]l5| , [Tr| , 

G n (r,A<)= (3) 
J dr'( [u(r' + r, t, At) - («>] [u(r', t, At) - («>] ), 

where. 



»(r, t, At) = IH(t, At) 5(r - r<(t)) . 
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G u (r,At) measures correlations in fluctuations of local 
displacements away from their average value. We are 
considering an equilibrium liquid and so G u does not de- 
pend on the choice of the reference time i. Similarly, 
(u) = (u(r,t,At)) does not depend on t; for a homoge- 
neous liquid, it also does not depend on r. In analogy to 
the relation between (n) and (N), we define the "total 
displacement" U(t,At) = J dru(r,t,At) and its ensem- 
ble average (U) = (U(t, At)). In a constant-iV ensemble, 
both (u) and (U) are readily evaluated from the mean 
displacement Jl = (iV _1 Yli=i M»(*> At)): (u) = ~p(n) and 
(U) = Jt(N). In equilibrium, (u), (U) and Jl do not de- 
pend on t, but they retain a dependence on At. 

G„(r, At) can be written so as to identify a spatial 
correlation function g u (r, At) analogous to g(r): 



correlation iimcuon g u \r, analogous 10 yyi). 

G„(r,At) = (N)J?6(r) + (u)(U) [<? u (r,At) - 1], 
where 
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The mean squared displacement p? = 
(iV^EiLi /•*!(*, At)), and also depends on At. The 
Fourier transform of G u (r, At) gives a "structure factor" 
S u (q, At) for the particle displacements: 5 u (q, At) = 

(^"^iliEjli ^(*.AtVi(*,Ai)exp[-iq. (r,(t) - 
In analogy to Eq. ||, the fluctuations of [/ are related 
to the volume integral of G u (r, At): 



[U- (C/)] 2 ) - y drG u (r,At) = (u)(U)kT Kv 
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We have defined the quantity n u in analogy to k. Hence, 
as for G(r), we can determine the large r behavior of 
G u (r, At) from the fluctuations of a bulk quantity, U. 
To evaluate these quantities we use data obtained 
from a molecular dynamics simulation of a model 
Lennard-Jones glass-former. The system is a three- 
dimensional binary mixture (80:20) of 8000 particles in- 
teracting via Lennard-Jones interaction parameters p7| . 
We analyze data from seven (p, P, T) state points on a 
line in the P, T plane approaching T c rj 0.435 at a pres- 
sure P w 3.03 |lq| . (In the remainder of this paper, all 
values are quoted in reduced units ]rj| ). The highest 
and lowest T state points simulated are (p = 1.09, P = 
0.50, T = 0.550) and (p = 1.19, P = 2.68, T = 0.451). 
Following equilibration at each state point, the particle 
trajectories are monitored in the NVE ensemble (E is 
the total energy) for up to 1.2 x 10 4 Lennard-Jones time 
units (25.4 ns in argon units) for the coldest T. Complete 
simulation details may be found in Q . All quantities pre- 
sented here are calculated using all 8000 particles in the 



liquid. The results presented here do not change when 
the minority particles are excluded @. 

For all seven state points, a "plateau" exists in both 
p? and the self part of the intermediate scattering func- 
tion F s (q,t) as a function of t M. The plateau separates 
an early time ballistic regime from a late time diffusive 
regime, and indicates "caging" of the particles typical of 
low T, high p liquids. The a-relaxation time r Q describes 
the decay of F s (q, t) to zero at the value of q correspond- 
ing to the first peak in the static structure factor S(q). 
Over the range of T studied, r Q increases by 2.4 orders 
of magnitude, and follows a power law r Q ~ (T — T c ) -7 , 
with T c = 0.435 and 7 ~ 2.8. The diffusion coefficient D 
follows a power law D ~ (T - T c )t, with T c = 0.435 and 
7 ~ 2.13, and thus diffusion and structural relaxation are 
"decoupled" |po| . The simulated liquid states analyzed 
here therefore exhibit the complex bulk relaxation behav- 
ior characteristic of a supercooled liquid approaching its 
glass transition. Both g(r) and S(q) for this liquid have 
been calculated previously |§,|l8|, and it has been shown 
that as T decreases, no long range structural correlations 
due to density or composition fluctuations occur. 




FIG. 1. g u {r, At) and g(r) versus r at T = 0.451. 
chosen on the order of r a . Inset: T(r, At) versus r. 
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In Fig.|l| we show g u (r, At) as a function of r for T = 
0.451, and with At chosen to be on the order of r a . g(r) 
for the same T is also shown. For a fixed choice of At, 
note that if the displacement were always the same for 
every particle, then g u (r, At) and g(r) would be identical 
for all r. Hence, it is deviations of g u (r,At) from g(r) 
that will inform us of spatial displacement correlations 
in excess of those that would be expected based on a 
knowledge of g(r) alone. We find that for this choice of 
At, g u (r,At) is appreciably higher than g(r) for values 
of r up to several interparticle distances. This excess 
correlation is made clearer in the inset of Fig^, where 
we show the function T(r, At) = [g u (r, At)/g(r)] — 1. 
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However, the question arises as to how to select the 
value of At. We find that the behavior of the liq- 
uid itself suggests a unique choice for At. To demon- 
strate this, we show in Fig. [2] the total excess correlation 
A = J drT(r,At) as a function of At. We find that 
there is a value of At = At* at which A is a maximum 
and that both the maximum value of A and At* increase 
with decreasing T. Hence for each T the spatial cor- 
relation of particle displacements is most prominent at 
At*. Moreover, all curves for T < 0.525 collapse onto 
a single master curve when t is scaled by At* and A is 
scaled by A(At*), suggesting that At* is a characteristic 
time for this liquid. In the remainder of this Letter, all 
quantities are therefore evaluated for At = At*. Fig. || 
shows that At* follows a power law with T: an excellent 
fit |l| to the form At* ~ (T - T c )~7 is obtained when 
T c = 0.435, and yields 7 = 2.3 ± 0.2. This value for 7 is 
different from the exponent found for r Q , but (within our 
numerical uncertainty) cannot be distinguished from the 
exponent governing the apparent vanishing of D at T c . 

If A is largest at At*, then we might also expect n u 
to be largest at At*, since by Eq. [?] k u quantifies the 
total magnitude (integrated over space) of the displace- 
ment correlations quantified by G u (t, At). We evaluate 
k u from the fluctuations of U according to Eq. fj] (Fig. ||b) 
and confirm n u exhibits the same behavior as A: k u goes 
to zero at short and long times, and has a maximum at 
a T-dependent At* . 




numerical error of T c , demonstrating that G u (r, At*) is 
becoming increasingly long-ranged as T — > T c . 
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FIG. 3. At* plotted versus T. The solid curve is a 
power-law fit to the data. INSET: Log-log plot of At* versus 
T — T c , and the power-law fit to the data. T c = 0.435. 
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FIG. 4. n u {At*) plotted versus T. The solid curve is a 
power-law fit to the data. INSET: Log- log plot of n u (At*) 
versus T — T c , and the power-law fit to the data. T c — 0.435. 



FIG. 2. (a) A versus At for different T. (b) k u as a function 
of At for the same T as in (a) . 

In Fig. |] we show the T-dependence of k u for At = 
At*. We find that n u (At*) grows monotonically with 
decreasing T , indicating that the range of the correlation 
measured by G u {r,At*) is growing with decreasing T. 
We find that a power law n u (At*) ~ (T - T c )" 7 fits well 
to the data when T c = 0.435, and gives 7 = 0.84. Thus 
k u exhibits an apparent divergence at a T that is within 



To estimate a correlation length associated with these 
displacement correlations, we evaluate S u (q,At*) for 
different T (Fig. . For intermediate and large q, 
S u (q,At*) coincides with S(q). However, for q — > 
a peak develops and grows with decreasing T, again 
demonstrating the presence of long range dynamical cor- 
relations. No growing peak at q — appears in the static 
structure factor S(q) (Fig. ||, inset). To quantify the cor- 
relation length, we attempted to fit S u (q, At*) using an 
Orstein-Zernike form, S u (q) oc 1/(1 + (, 2 q 2 ), where £ is 
the correlation length. Although this form fits well to 
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the data at the highest T, it fails at lower T, making the 
interpretation of the fitted £ values ambiguous. 

Nevertheless, it was shown previously for this system 
that highly "mobile" particles move cooperatively and 
form clusters Q whose mean size diverges at T c |8| . These 
clusters contribute to the growing range of G u (r, At*), 
and thus they can be used to give a rough estimate of the 
length scale over which particle motions are correlated. 
At T — 0.451, this average length scale exceeds 3 particle 
diameters, and the largest cluster has a length scale that 
exceeds the size of our simulation box (approximately 19 
particle diameters on a side.) 
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FIG. 5. S u (q, At*) versus q for different T. The values at 
q — are obtained from the fluctuations in U via the relation 
S u (q = 0, At*) = pTn u (At*). INSET: Static structure factor 
S(q) for four different T. 

In summary, we have defined a correlation function 
that quantifies the spatial correlation of single-particle 
displacements in a liquid. Using this function, we have 
shown in computer simulations of an equilibrium liquid 
that the displacements of particles are spatially corre- 
lated over a range and time scale that both grow with 
decreasing T as the mode coupling temperature is ap- 
proached. While MCT makes no predictions concerning a 
growing dynamical correlation length [^,^3) , calculation 
of the vector displacement-displacement correlation func- 
tion may be tractable within the mode coupling frame- 
work. We have also identified a bulk dynamical variable 
U whose fluctuations appear to diverge at T c . Hence, U 
is behaving much like a static order parameter on ap- 
proaching a second-order phase transition. Our analy- 
sis therefore suggests that an extension to dynamically- 
defined quantities of the framework of ordinary critical 
phenomena may be useful for understanding the nature 
of supercooled, glass-forming liquiods. 
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